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Using the Bethe ansatz we obtain the exact solution for the one- 
dimensional asymmetric avalanche process. We evaluate the velocity of dis- 
persive flow as a function of driving force and the density of particles. The 
obtained solution shows a dynamical transition from intermittent to continu- 
ous flow. 
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I. INTRODUCTION AND THE MODEL DEFINITION. 

The avalanche dynamics is a basic scenario of relaxation of unstable states in extremal 
systems where each movable element is near a border of stability. A typical long-tailed 
distribution of avalanche sizes leads to the dispersive transport of particles KJ. As an 
illustrative example, granular systems exhibit intermittent avalanches which enables one 
to use granular piles (sand piles, rice piles) for explanation of self-organized criticality in 
generic dissipative systems ||. In the past decade, it has become clear that the dispersive 
transport can be recast in terms of interface depinning []3],|4j and various growth models || . 
Recently, a dynamical transition from intermittent to continuous flow in a random sandpile 
model has been revealed 0. Nevertheless, an explicit theoretical description of stochastic 
avalanche processes and an exact evaluation of characteristics of dispersive flow remains an 
open problem. 

Despite drastic simplifications which were introduced to mimic real avalanches, exact 
results are scarce even for the determimistic dynamics. As to stochastic dynamics, it is es- 
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pecially difficult as it is beyond the class of abelian models |7j, where asymmetric processes 
appear to be solvable ||. The situation is to be compared with the theory of exclusion pro- 
cesses where many properties, such as steady states, average current, diffusion constant etc. 
have been calculated for an asymmetric one-dimensional case p|^T2| . The usual presentation 
of the asymmetric exclusion process (ASEP) is given by a master equation for the probabil- 
ity Pt(xi, ...,xp) of finding P particles at time t on sites xi, ...,xp of a ring consisting of N 
sites. During any time interval dt, each particle jumps with probability dt to its right if the 
target site is empty. This elementary restriction leads to a non-trivial problem of evaluation 
of the steady state properties, which can be solved by the Bethe ansatz. 

In a similar way, the simplest asymmetric avalanche process (ASAP) can be formulated 
as follows. In a stable state, each of iV sites on a ring is either occupied by one particle 
or empty. The total number of particles P is fixed. During time interval dt, each particle 
jumps with probability dt to its right. In the course of time, some site x may get unstable 
with occupation number n > 1. Then it must relax immediately to the stable state by 
transferring to its right either n particles with probability fi n or n — 1 particles with the 
probability 1 — /i n . The quantity \i n can be associated with a driving force acting on the 
unstable group of n particles. 

The main difference between the ASEP and ASAP lies in the depth of reconstruction of 
a configuration C = {x\, xp} during the time interval dt. In the ASEP, the total distance 
Y t covered by all particles between time and t increases by 1 during dt if the configuration 
C differs from a new one C or remains unchanged if the motion is forbidden. In the ASAP, 
the motion of a particle is always possible and increase of Y t is not bounded. Thus, the 
configuration C may be completely different from C' depending on numbers of particles 
spilled to right from each unstable site. 



The present formulation of the ASAP is inspired by works where a model of activated 



random walks is introduced and [[DJ where the directed avalanche dynamics is formulated 
in terms of continuous variables. Under an assumption about independence of variations 
of the avalanche size at each time step, the probability distribution of avalanche sizes is 



found exactly fL5| , [0. However, the configurational space of the continuous model is too 
complicated to determine steady state features. 

Here, extending the Bethe anzatz approach to exclusion processes, we obtain the expres- 
sion for the generating function of Y t for the discrete ASAP in the thermodynamic limit of 
large N for a fixed density of particles p = P/N. We find two phases corresponding to a 
dispersive flow and a continuous flow, and evaluate the exact average velocity in the whole 
range of parameters of the first phase. We determine the separation line between two phases 
where avalanches are critical. 



II. THE DYNAMICAL RULES AND THE AVERAGE VELOCITY. 

Even before going into details of our calculations, we can get important restriction on 
the toppling rules of the ASAP. For the problem to be solvable by Bethe ansatz, one has 
to make sure that many-particle problem can be reduced to the problem of two interacting 
particles. In other words, the toppling probabilities /i n of unstable configuration with n > 1 
particles at the same site should be determined recursively in terms of two-particle toppling 
probability /i2 only. Namely, the probability for n particle to leave the site is the sum of two 
processes. In the first process two particles leave the site with probability fi2 and then n — 2 
particles leave the site with probability p n -2- The second process corresponds to spilling two 
particles with probability 1 — /i2 when only one particle leave the site and then remaining 
n — 1 particles leave the site with probability p n -i- Thus, we obtain the recursion relations 
which express all the probabilities trough the only constant /i 

fJ.i = 
A*2 = fi, 

Hn= Hn-2H + Hn-\(1 — H), (1) 

or in the form of one step recursion, 

/i n = - /i n _i) (2) 
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Although the exact solution of the Bethe anzats requires a long technical analysis, the 
average velocity of the dispersive flow can be obtained with minimal assumptions from 
simple combinatorial arguments. The average velocity of the particle flow in the ASAP is 
determined by the average number of steps of all particles involved into an avalanche during 
the time interval t and can be written as 

The only assumption we use is that the probability for any site of the infinitely large lattice 
to be occupied does not depend on state of the other sites and is equal p. Let us consider 
the avalanche starting at the site i. To calculate the velocity, we introduce the probability 
Pi,j(n) for n particles to be transferred from the site j to its right provided the avalanche is 
started at site i. Defining also the total probability to transfer exactly n particles from any 
site during the whole avalanche P{n) = Y^jLiPi,j{n) we can express the average velocity as 
follows 



oo 

r 

71=1 



= £ "P(n) (4) 



Using the translation invariance of the stationary state one can rewrite P{n) as the sum of 
topplings of one site % in all avalanches started on the left of the site %. 

oo 

P(n) = J2^An). (5) 
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The dynamic rules of the model relate the values of P(n) for two neighboring sites. Due to 
the translation invariance, P(n) does not depend on site, and obeys the recurrent relations 

P(n) = P{n - l)pp n + P(n)(p(l - p n+l ) + (1 - p)p n ) + P{n + 1)(1 - p(l - p n+1 )). (6) 

where we express P(n) at the site i via P(n) at the site % — 1. For P(l) we have different 
recurrent relation. The term with P(n — 1) should be replaced by 1 which corresponds to 
the case when avalanche is starting at given site, 

P{1) = 1 + P(l)(p(l - p 2 )) + P(2)(l - p)(l - p 2 ). (7) 



To find the solution of these recurrent relations one should also fix the value of P(l). If we 
want the system to be stationary we should make sure that the number of starting avalanches 
is equal to the number of dying ones. There is P(l)(l — p) avalanches ending at any site for 
each avalanche starting at the same site, what gives P(l)(l — p) = 1. Then the solution of 
the system of the relations ([], 0) is 

P(n) = / /" P(n- 1), (9) 
1 - p 1 - /i n 



or applying these relations recursively 
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P\ 1 -PJ t=2 l ~ Pr. 

Substituting the expression for P{n) into Eq.(^), and using the recursion (g), we get 
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^-t — iT*-) (11) 

As we shall see below, this formula coincides with that we obtain from the exact Bethe 
anzats solution. This means that our assumption about uncorrelated stationary state is 
valid in the thermodynamic limit. 

III. BETHE ANZATS SOLUTION. 

Consider the ASAP consisting of P particles on a ring of sites and denote by Pt{C) 
the probability of finding at time t the system in a configuration C. The probability Pt(C) 
satisfies 

^Pt(C) = J2[M (C,C') + M 1 (C,C']P t (C') (12) 
at c , 

where Mi(C, C')dt is the probability of going from C to C during the time interval dt, and 
Mq is a diagonal matrix 

M (C,C) = - £ Mi [d , C) (13) 

c'+c 



Before using the Bethe ansatz, it is instructive to note that in the region where the distances 
between every two neighbouring particles exceed 1, the master equation ( |T2"D becomes "free": 

■^-P t (xi, x P ) = 53 e 7 Pt(xi, x k - 1, ...,x P ) - P t (xi, x P ) (14) 

where exp(7) is activity of a single step. To compensate the difference between (Tl2|) and 



(|14]) when — x^-i = 1 for some fc, we introduce the boundary conditions 

P t {...,x,x,...) =e^{l-fx)P t {...,x-l,x,...) + e 2 ^P t {...,x-l,x-l,...) (15) 

This condition can be viewed as the recurrent relation where the "intermediate" probability 
of an unstable configuration Pj(..., x, x, ...) is given in terms of another unstable configuration 
P t (...,x — l,x — 1, ...) and so on. All boundary conditions for more then two particles can 
be reduced to the two-particle case. This implies a recurrent relation for the probability fi n 
which is nothing but the two particle reducibility (@) discussed above. 

Now, we can define the ASAP by ( |i~4"D and fll5|) instead of fli~2|) without even knowing 
the exact form of the matrix M(C, C') which is very cumbersome for the ASAP model. 
Specifying a configuration C by positions 1 < x± < x^--- < xp < N of the P particles, we 
use the Bethe ansatz for an eigenvector of the matrix M + Mi in the form 

E^n^a as) 



" °QU) 
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where the sum is over all of the permutations Q of 1, 2, P. The condition ( pT5[) fixes the 
two particle S-matrix A^/A,-, as 

A jk _ 1 - (1 - Li)eiZj - fie^ZjZk 



A kj 1 - (1 - fi)e"fz k - /le^ZjZk 
Imposing the periodic boundary conditions gives the Bethe equations 



(17) 
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-N 



h 



— 1 — (1 — ju)e 7 2j — fxe 21 'zjZk 



H)e<z k - fxe^ZjZk 



(18) 



The eigenvalue A (7) corresponding to (fTB) is 



A( 7 ) = -P + e 7 ^^ (19) 



The dependence of the eigenvalue on 7 allows one to use A (7) as the large deviation function 
of Y t . Specifically, the average velocity of the particle in ASAP is expressed through the 
derivative of A (7). 

1 dA(7) 



P dq 



(20) 
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IV. THE BETHE EQUATIONS IN THERMODYNAMIC LIMIT. 

The Bethe ansatz equations (|Tj|) together with ([19]) give the exact solution of the problem 
for all N and P. The rest of the paper is devoted to evaluation of v in the thermodynamic 
limit N — > 00 for a fixed density of particles p = P/N. For a finite N, the largest eigenvalue 
A corresponds to the solution {zj} which converges to Zj = l,j = 1, ...,P as 7 — > 0. For 
small 7 > 0, the distance \zj — 1| grows rapidly with N for all j and becomes of order of 1 
in the limit N — > 00. Introducing a variable a by 

1 _ P ia j 

Zj = - e~i (21) 

and assuming the solutions {<yj} are distributed along a smooth curve in the complex plane 
a = (u + ir) with endpoints (—a + ib) and (a + ifo), we obtain the Bethe equation in the 
form 

p(a) = 27rF(a) + — ( a+lh 6 (a - P)R(p)df3 - z 7 (22) 

2ll J-a+ib 

where we defined as usual a function F(a) such that dF/da = —R(a)/2ir and F(— a + ib) = 
—F(a + ib) = p/2. The functions p(a) and 8(a) are 

^) = - iln (lT^) (23) 

and 

a/ \ -i /Cosh(z/ + ia/2) 

0(a) = -2 In — 24 

cosh(z/ — za/2) 

where z/ = — ln(/i)/2. Taking the derivative in (p2|), we get the integral equation for i?(a) 



with 



and 



R(u, b) + — f a K(u - v)R(v, b)dv = £(u, b) (25) 

Z7T J -a 



j., \ coshz/ . . 

sinn v — smn(v — 



r w x sinh2i/ , . 

AT(a) = — 27 

cosh 2u + cos a 



All that is very similar to the equations for the asymmetric 6- vertex model [ 17 , |18| (see 



also fL9| ) with an essential exception: both terms containing Zj and ZiZj in (18) are negative, 
which is the reason for a dynamical transition, as we shall show below. 



If a = 7T, equation (|25[) can be solved by the Fourier transformation. To evaluate <9 7 A, 
we have to find the solution of (|25| ) in a vicinity of the point a = tt which corresponds to 
a "conical" point, considered in |Tj|. Following Bukman and Shore, we write the solution 
R(u) as an expansion in e = n — a up to order of 0(e 4 ) 

R(u) = Ro(u) + e^Rxiu) + e 2 5R 2 (u) + e 3 8R 3 (u)e*5R 4 (u) (28) 

The necessity of such a long expansion will be seen in further calculations. The Fourier 
transformation is defined by 

oo 

X(u) = Yl (X) n e~ mu (29) 

n=— oo 

where X stands for Rq, 5R m , £, K. The non-zero Fourier coefficients of K and £ for b > —1v 



arc 



{K) n = (_i)« e - 2 " |n| (30) 
(O n = -e k (l-(-l)V w ),n<0 (31) 



Then, (|25| ) gives 
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1-e 

In the next order in e, we have 



R (u) = (i2o)o + - — - b (32) 



(8R 1 ) n (K n - 1) = ^°M(-l)»(jr)„ (33) 

7T 



arc 



so, that Ro(tt) = 0,(8Rx) n = and (Rq)o = 1/(1 + expfr). The next terms in (|2"8" 
evaluated in |Tj| 

^N = -k'W (34) 
o 



,«,(„) = _ E MM^ e -*- (35) 
n ^ 3tt(1 - [K) n ) 

8R A = -^Rf\<K) (36) 

Thus, in the expansion of R(u), 5R 2 (u) = 5R 2 and 8R±(u) = 5i? 4 are real constants and 
5R^{u) is imaginary. 



Now, we are ready to start a direct evaluation of <9 7 A in (plj) using <9 7 A = <9 a A/<9 a 7. First 



we find da'y. To this end, we put a = a + ib in (|22|), and take the derivative by a at a = it. 
Recalling the conditions F(a + ib) = —p/2 and #(0) = 0, we obtain 

ird a p + z<9 7 = R(a, b) + (2n)- 1 9(2a)R(-a, b) + 

(27T)- 1 J* a 6{a-v)d a R{v,b)dv (37) 

To evaluate r.h.s. of (|37j), we express the values R(a,b),R(—a,b) and #(2a) by their Tailor 
expansions at a = ir up to the order of 0(e 3 ) using fl2~j|) and (^). The integral in (^) is 
treated as 

/(v)dt;= f f{v)dv + B(e) (38) 

J —Tx+e J —IT 



with 

oo 

M! 



m=l 



and, therefore, can be evaluated by the Fourier transformation. The Fourier coefficients of 
0(n — v) are 

27T? 

(0(tt - v)) n = — (e~ 2 ^ - (-l) n ),n ^ (39) 

7?; 



and (9(tt — v)) = 2ir 2 . The only f-dependent part of d a R(v,b) is 8R 3 (v). Using ( |3"0"| , |3"4| - |3T)| ) . 
and (|39f) , we obtain the explicit expression for r.h.s. of ( |3~7|) 

^W-^ + ^flfW (40) 

Due to (|32|) , the first and third terms in r.h.s. of ( |40"D are real, the second one is imaginary. 
Therefore, we have 

d al = *e 2 ^ (41) 

7T 

The expression for d a X can be found in a similar way. In this case, we take the derivative 
by a in 

± = -L [ a R(u,b)(z(u,b)-l)du (42) 
where z(u, b), according to Q2I|), * s 



*M) = TT^=S (43) 



j ^iu—b 

1 + e* 

The obtained derivative is similar to ( |37D but contains z(v, 6) instead of 6*(a — v). So, we 
need the Fourier coefficients of z(v, b) which are 

(z) n = (1 + e 2u )(-l) n e {2v+h)r \ n < (44) 



and (2)0 = 1. Continuing as in (|37|), we get 

1 <9A 2 i?o(7r)z'(7r, 6) 



V < )( > 7T 



3e 2 £(z)_ n («? 3 )n (45) 



We can see that the term 5R^{u) in is relevant. As to SR^u), it is sufficiently that it 
is a constant and does not lead to a divergency by integration. 
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Substituting the explicit expressions for (z)- n and (5Rs) n gives the second term in r.h.s. 
of fl45| ) in the form 

7'/?' (V) 00 (—]) n p-( 4u + b ) n 

^ {1 + e ^3M^UlJ__ ,46) 

Due to (H), i?o (7r) is cancelling in d a \/d a ^. Then, using and the identity p = (-Ro)o — 
1/(1 + exp b) we obtain the final result 

_ (i- P )(i +/i ) i+p - (-lrv" / P r 
(i-p(i+^)) 2 w ^i-(-i)"^ u-py ' 

which exactly coincides with the formula (^T|) after some algebra. 

The velocity of flow v diverges at p c = 1/(1 + p) which implies a transition to the phase of 
continuous flow. The value of critical density p c can be easily understood from the condition 
of a balance between gaining (p c ) and losing (1 — q^) one particle at each step of a large 
avalanche. 

The considered model is a directed version of the model of activated random walks intro- 
duced in |li| to see how a conservative dynamical system with the sandpile toppling rules 
approaches criticality. It has been shown in [14| that the relaxation time r and correlation 
length £ diverge as r ~ |p c — p\~ Ul and £ ~ |p c — p\~ U2 . The exponents v\ and z/ 2 have been 
determined numerically for several kinds of toppling rules. In the directed case, £ coinsides 
with r and is proportional to the average size of avalanches (s). On the other hand, (s) = v 
so we have from PT| ) (s) ~ (p c — p) -2 and Ui = u 2 = 2. An extension of this result to the 
symmetrical case []2D| is a very interesting and difficult problem. 
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